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Abstract. 

The instantaneous state of a neural network consists of both the degree of excitation 
of each neuron, the network is composed of, and positions of impulses in communication 
lines between neurons. In neurophysiological experiments, the neuronal firing moments 
are registered, but not the state of communication lines. But future spiking moments 
depend essentially on the past positions of impulses in the lines. This suggests, that the 
sequence of intervals between firing moments (interspike intervals, ISIs) in the network 
could be non-Markovian. 

In this paper, we address this question for a simplest possible neural "net" , namely, 
a single neuron with delayed feedback. The neuron receives excitatory input both from 
the driving Poisson stream and from its own output through the feedback line. We 
obtain analytical expressions for conditional probability density P(t n _|_i | f„, . . . , ti, to), 
which gives the probability to get an output ISI of duration t n+ \ provided the 
previous (n + 1) output ISIs had durations t n , . . . , t±, to- It is proven exactly, that 
P(i„+i | t n , . . . , ti, to) does not reduce to P(t n+ i \ t n ,...,tx) for any n > 0. This 
means that the output ISIs stream cannot be represented as Markov chain of any 
finite order. 

PACS numbers: 87.19.11, 87.10.-e, 02.50.Cw, 02.50.Ey, 87.10. Ca, 87.10.Mn 
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1. Introduction 

In a biological network, the main component parts are neurons and interneuronal 
communication lines - axons p]. These same units are the main ones in most types of 
artificial neural networks [2]. If so, then the instantaneous dynamical state of a network 
must include dynamical states of all neurons and communication lines the network is 
composed of. The state of a neuron can be described as its degree of excitation. The 
state of a line consists of information of whether the line is empty or conducts an impulse. 
If it does conduct, then further information about how much time is required for the 
impulse to reach the end of the line (time to live) describes the line's state. 
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Figure 1. Signal processing in the binding neuron model [4]. 

In neurophysiological experiments, the triggering (spiking, firing) moments of 
individual neurons are registered. The sequence of intervals between the consecutive 
moments (interspike intervals, ISIs) is frequently considered as renewal stochastic 
process. Recently, based on experimental data it was offered that the ISIs sequence 
could be Markovian of order 4 or higher [3]. 

The presence of memory in the ISI sequence is not surprising, taking into account 
that information about triggering moments leaves unknown the states of communication 
lines at those moments. On the other hand, it is namely the impulses propagating 
in the communication lines that connect past firing moments with the future ones 
in a reverberating neural network. Without knowledge of communication line states, 
information about previous neuronal firing moments could improve our predicting ability 
of the next ones. The exact answer of what kind of memory could be expected in an 
ISI sequence of a neuron embedded in a reverberating neural network driven with some 
noisy stimulation requires rigorous mathematical treatment. 

In this paper, we consider a simplest neural "net", namely, a single neuron with 
delayed feedback, which is driven with Poisson process. As neuronal model we take 
binding neuron as it allows rigorous mathematical treatment. We study the ISI output 
stream of this system and prove that it cannot be presented as Markovian chain of any 
finite order. This suggests that activity of a more elaborate network, if presented in 
terms of neuronal firing moments, should be non-Markovian as well. 

2. The object under consideration 

2.1. Binding neuron model 

The understanding of mechanisms of higher brain functions expects a continuous 
reduction from higher activities to lower ones, eventually, to activities in individual 
neurons, expressed in terms of membrane potentials and ionic currents. While this 
approach is correct scientifically and desirable for applications, the complete range of 
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Figure 2. Binding neuron with feedback line under Poisson stimulation. Multiple 
input lines with Poisson streams are joined into a single one here. A is the delay 
duration in the feedback line. 

the reduction is unavailable to a single researcher /engineer due to human brain limited 
capacity. In this connection, it would be helpful to abstract from the rules by which 
a neuron changes its membrane potentials to rules by which the input impulse signals 
are processed in the neuron. The coincidence detector, and temporal integrator are the 
examples of such an abstraction, see discussion in [5]. 

One more abstraction, the binding neuron (BN) model, is proposed as signal 
processing unit |6], which can operate either as coincidence detector, or temporal 
integrator, depending on quantitative characteristics of stimulation applied. This 
conforms with behavior of real neurons, see, e.g. [TJ. The BN model describes 
functioning of a neuron in terms of discret events, which are input and output 
impulses, and degree of temporal coherence between the input events, see Figure [TJ 
Mathematically, this is realized as follows. We expect that all input impulses in all 
input lines are identical. Each input impulse is stored in the BN for a fixed time, r. 
The t is similar to the tolerance interval discussed in [Sj. All input lines are excitatory. 
The neuron fires an output impulse if the number of stored impulses, S, is equal or 
higher than threshold value, N . After that, BN clears its memory and is ready to 
receive fresh inputs. That is, every input impulse either disappears contributing to a 
triggering event, or is lost after spending r units of time in the neuron's internal memory. 
It is clear, that BN fires when a bunch of input impulses is received in a narrow temporal 
interval. In this case the bunch could be considered as compound event, and the output 
impulse — as an abstract representation of this compound event. One could treat this 
mechanism as binding of individual input events into a single output event, provided 
the input events are coherent in time. Such interpretation is suggested by binding of 
features/events in largescale neuronal circuits [9| [T0| [TTj. 

Further, we expect that input stream in each input line is the Poisson one with some 
intensity Aj. In this case, all input lines can be collapsed into a single one delivering 
Poisson stream of intensity A = J2i see Figure El 

For analytical derivation, we use BN with N = 2. The case of higher threshold is 
considered numerically. 
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2.2. Feedback line action 

In real neuronal systems, a neuron can have synaptic connection of its axonal branch 
at its own dendritic tree, see [121 [13] f° r experimental evidence. As a result, the neuron 
stimulates itself obtaining excitatory impulse after each firing with some propagation 
delay. We model this situation assuming that output impulses of BN are fed back 
into BN's input with delay A. This gives BN with delayed feedback, Figure [2j See 
also Supplementary Matherial for animation of BN with delayed feedback in action. 
Impulses from the feedback line have the same excitatory action on BN as those arrived 
from Poisson stream. Namely, each one of them is stored in BN's memory for time r, 
after which it desappears completely, see section 12.11 

The feedback line either keeps one impulse, or keeps no impulses and cannot convey 
two or more impulses at the same time. If the feedback line is empty at the moment of 
firing, the output impulse enters the line, and after time interval equal A reaches the 
BN's input. If the line already keeps one impulse at the moment of firing, the just fired 
impulse ignores the line. 

Any output impulse of BN with feedback line may be produced either with impulse 
from the line involved, or not. We assume that, just after firing and sending output 
impulse, the line is never empty. This assumption is selfevident for output impulses 
produced without impulse from the line, or if the impulse from the line was involved, 
but entered empty neuron. In the letter case, the second (triggering) impulse comes 
from the Poisson stream, neuron fires and output impulse goes out as well as enters 
the empty line. On the other hand, if impulse from the line triggers BN, which already 
keeps one impulse from the input stream, it may be questionable if the output impulse 
is able to enter the line, which was just filled with another impulse. We expect it does. 
This means that the refraction time of biological axon modelled as feedback line is equal 
A. Thus, at the beginning of any output ISI, the line keeps impulse with time to live s, 
where s e]0; A]. In this paper, we consider the case 

A<r (1) 

in order to keep expressions shorter. 

3. Statement of the problem 

The input stream of impulses, which drives neuronal activity is stochastic. Therefore, 
the output activity of our system requires probabilistic description in spite of the fact 
that both the BN and the feedback line action mechanisms are deterministic. We treat 
the output stream of BN with delayed feedback as the stationary process^!- In order to 
discribe its statistics, we introduce the following basic functions: 

| The stationarity of the output stream results both from the stationarity of the input one and from 
the absence of adaptation in the BN model, see Section 12.11 In order to ensure stationarity, we also 
expect that system is considered after initial period sufficient to forget the initial conditions. 
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• joint probability density P(t m , i m _i, . . . , to) f° r ( m + 1) successive output ISI 
durations. 

• conditional probability density P(t m | i m _i,...,t ) for output ISI durations; 
P(t m | t m _i, . . . , to)dt m gives the probability to obtain an output ISI of duration 
between t m and t m + dt m provided previous m ISIs had durations t m _i, t m _ 2 , • • • , to, 
respectively. 

Definition T/ie sequence of random variables {tj}, taking values in Q, is called the 
Markov chain of the order n > 0, if 

V m >nVt og Q • • • ^t m efl P\pm I tm— lj • • • , ^0) = P\pm I t m _i, . . . , t m _ n ), 

and t/w's equation does not hold for any n' < n (e.g. JXffiJ. In the case of ISIs one 
reads Q = R + . 

In particular, taking m = n + 1, we have the necessary condition 

P(t n+ 1 | i n , • • • MM) = P(tn+1 \t n , ■ ■ ■ M), U G O, Z = 0, . . . , 71 + 1, (2) 

required for the stochastic process {tj} to be the n-order Markov chain. 

Theorem 1 The output ISIs stream of BN with delayed feedback under Poisson 
stimulation cannot be represented as a Markov chain of any finite order. 

4. Proof outline 

In order to prove the Theorem 1, we are going to show analytically, that the equality (J2J) 
does not hold for any finite value of n, namely, in the exact expression for conditional 
probability density P(i n +i | t n , . . . ,ti,t ), elimination of i -dependence is impossible. 
For this purpose we introduce the stream of events (t, s) 

ts = {. . . , (tj, Si), . . .}, 

where S{ is the time to live of the impulse in the feedback line at the moment, when 
ISI tj starts. We consider the joint probability density P(t n+1 , s n+ i; t n , s n ; . . . ; to, «o) 
for realization of (n + 2) successive events (t, s), and the corresponding conditional 
probability density P(t n+ i, s n+ i \ t n , s n ; . . . ; t , sq) for these events. 

Lemma 1 Stream ts is 1-st order markovian: 
V„>oV to> oV soe ] 0; A] • • • Vi n+1 >oV Sn+ie ] ;A] 

P(in+1, s n+l | t n , S n ] . . . ; to, So) = P(t n +i, S n +i I t n , s n ), 

where {t , . . . , t n +i} is the set of successive ISIs, and {s , 
corresponding times to live. 

Proof Indeed, the value of s n+ i characterizes the state of the system at the moment of 
triggering, 6, and the value of t n+ i characterizes the system's behavior after that 
triggering, which means that, in physical time, s n+ i always gets its value before 
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than the t n+ i does. Once the value of s n+ i is known, the realization of t n+ \ is 
completely determined by a unique realization of the input Poisson process after 
the 9. 

At the same time, in P(t n+1 , s n+ i \ t n , s n , . . . , t , s ) the value of s n+ i can be derived 
unambiguously from (t n ,s n ) (See Sections 12.21 and l5\2l : 

Sn+l &n t ni t n <C S n , 

= A, t n > s n . (4) 

Just after triggering, BN appears in the standard state (it is empty), the state of 
line is given by the value of s n+ i, and the state of input Poisson stream is always 
the same. Therefore, once the pair of values (t n , s n ) is given, the state of the system 
at the moment of (n+ l)-th ISI beginning is determined completely, and knowledge 
of previous values of (U, Sj), i < n adds nothing to our predictive ability as regards 
the values of (t„+i, s n+ i), which proves 

In order to find the conditional probability density P(t n+1 \ t n , . . . , ti, t ), the 
following steps should be performed: 

• Step 1. Use property (J3j) for calculating joint probability of events (t,s): 

P\pn+li s n+l'i tni s n\ ■ ■ ■ i ^0) s o) = 

P{tn+lj s n+l | t n , S n ) . . . P(t h Si | t , S )P(t , S ), (5) 

where P(t,s) and P(t n ,s n \ t n _i,s n _i) denote the stationary probability density 
and conditional probability density (transition probability) for events (t, s). 

• Step 2. Represent the joint probability density for successive output ISI durations 
as marginal probability by integration over variables s^, i — 0, 1, . . . , n + 1: 

P(t n +i, t n , ■ ■ ■ , to) = 

r-A pA pA 

ds dsi... ds n+1 P(t n+ i,s n+ i;t n ,s n ; . . . ;t ,s ). (6) 



'o jo jo 
• Step 3. Use the definition of conditional probability density: 

P{t n +l | tn, ■ ■ ■ , h,t ) = — ' — ' — . (7) 

f[V n , . . . , to J 

Taking into account Steps 1 and 2, for joint probability density P(t n +i, ■ ■ ■ ,t ) one 
derives 

P(t n +i, t n , ■ ■ ■ , to) = 

l>A pA w+1 

/ ds ... ds n+ iP(to,s ) J\P(t k ,s k \t k -i } s k -.i). (8) 
Jo Jo k=1 

In the next section, we are going to find the exact analytical expressions for 
probability densities P(t,s) and P(t k ,s k \ tfc-i, Sfc-i), and perform the integration in 
OH]). Then we aply the Step 3, above, to find expressions for conditional probabilities 
P(t n +i | t n , . . . ,ti,t ). It appears, that the conditional probabilities have singular parts 
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Figure 3. Output ISI probability density P(t) (a) and probability density f(s) for 
times to live of the impulse in the feedback line (b), found analytically in [15]. Here r 
= 10 ms, A — 8 ms, A = 150 s _1 , Nq=2. The presence of 5-function in both densities 
is clearly visible. 



of the Dirac's 5-function type. This is because the system's dynamics involves discret 
events of obtaining impulse by neuron (see below). In order to prove that the equality (J2]) 
does not hold for any n > 0, we use the singular parts only. 

5. Main calculations 

5.1. Probability density P(t, s) for events (t, s) 

The probability density P(t, s) can be derived as the product 

P(t,s)=F(t\s)f(s), (9) 

where f(s) denotes the stationary probability density for time to live of the impulse in 
the feedback line at the moment of an output ISI beginning, F(t \ s) denotes conditional 
probability density for ISI duration provided the time to live of the impulse in the 
feedback line equals s at the moment of this ISI beginning. Exact expressions for both 
f(s) and F(t \ s) are given in [T5| Eqs.(5),(6) and (31)]. In this paper we need only 
singular parts of those expressions, which read: 

F sing (t | s) = \se~ Xs 5(t - s), (10) 

4p2AA 

/ smg (s) — a - 5(s — A), where a=- — r , (11) 

J K ' K h (3 + 2AA)e 2AA + 1' v ; 

where a gives the probability to obtain the impulse in the feedback line with time to live 
equal A at the beginning of an arbitrary ISI, A — is the input Poisson stream intensity. 

The presence of 5-functions in F(t \ s) and f(s) can be explained as follows. The 
probability to obtain an output ISI of duration t exactly equal s is not infmitesimally 
small. Due to ([[]), it equals to the probability to obtain exactly one impulse from the 
Poisson stream during time interval ]0;s[, which is Ase _As . The second impulse comes 
from the line and triggers the neuron exactly after time interval s. So, we have the 
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non-zero probability to obtain an output ISI of duration exactly equal to s. This gives 
the ^-function at t = s in the probability density F(t \ s). 

The probability to have time to live, s, exactly equal A at the moment of an output 
ISI beginning is not infinitessimally small as well. Every time, when the line is free at 
the moment of an output ISI beginning, the impulse enters the line and has time to 
live equal A. For the line to be free from impulses at the moment of triggering, it is 
nessesary that t > s for the previous ISI. The set of realizations of the input Poisson 
process, each realization results in t > s, has non-zero probability a, see ( fTTl) . and this 
gives the ^-function at s = A in the probability density f{s). 

The output ISI probability density P(t) can be obtained as the result of integration 
of © (see |I5] for details): 

P(t) = f F(t\s)f(s)ds. (12) 
Jo 

Examples of P(t) and f(s) graphs are given in Figure |3j 
5.2. Conditional probability density P(t k , Sk \ tk-ii Sfc-i) 

Here we find the conditional probability density P(t k , s k | t k _i, Sk-i) for events (t k , s k ), 
which determines the probability to obtain the event (t k ,s k ), with precision dt k ds k , 
provided the previous event was (tfc-i, Sk~i). By definition of conditional probabilities, 
the probability density wanted can be represented as the following product 

P(tk,Sk | tk-i,Sk-i) = F{t k | Sk,tk-i, Sk-i)f(sk | tfc_i,Sfc_i), (13) 

where F(t k \ s k , t k -i, s/c-i) denotes conditional probability density for ISI duration, t k , 
provided i) this ISI started with lifetime of impulse in the feedback line equal to s k , and 
ii) previous (t, s)-event was (tfc_i, Sfc-i); f($k \ tfc-i, Sfc-i) denotes conditional probability 
density for times to live of the impulse in the feedback line under condition ii). It is 
obvious, that 

F(t k | Si,t w ,vi) = F(t k | s k ), (14) 

because with s k being known, the previous event (t k _i,Sk-i) does not add any 
information, useful to predict t k (compare with proof of Lemma 1). 

In order to find the probability density f(s k \ t k _i, s k ^i), let us consider different 
relations between and s k -i- If tfc-i > Sfc-i, the line will have time to get free from 
the impulse during the ISI t k -i- That is why at the beginning of ISI t k , an output spike 
will enter the line and will have time to live equal s k = A with probability 1. Therefore, 
the probability density contains the corresponding delta-function: 

f(sk I *k-i,Sfc-i) = H s k - A )> > s fe -L. (15) 

If < Sfc-i; than the ISI t k _i ends before the impulse leaves the feedback line. 
Therefore, at the beginning of the t k , the line still keeps the same impulse as at 



Delayed feedback causes non-Markovian behavior of neuronal firing statistics 



9 



the beginning of tk-i- This impulse has time to live being accurately equal to 

s k = S fc _i - t k -l, so 

f(sk I tk-i, s fc _i) = 8(s k - s fc _i + t fc _i < s fe _i. (16) 

Taking all together, for the conditional probability density P(tk,Sk | ijt-i, sjt-i) one 
obtains 

P(t fc , s fc | tfc-i, Sfe_i) = | s k )5(s k - A), t fc _i > s k -i, 

= F(t k | s k )5(s k - s fc _i + t fc _i), t fc _i < s fe _i, (17) 
where exact expression for F(t | s) is given in [151 Eqs.(5), 



5.3. Joint probability density P(t n+ i, . . . , to) 

In this section, we are going to find the exact analytical expression for the joint 
probability density P(t n +i, ■ ■ ■ , to) at the domain 

D 1 = f(t ,...,t n ) | J><aJ. (18) 

It is worth to notice, that the set of (n + 1) successive ISI durations t ,...,t n has 
non-zero probability, p& > 0, to fall into the domain (j!8p . Indeed, BN with threshold 
iVo = 2 needs 2(n + 1) input impulses within time window ]0; A[ to be triggered (n + 1) 
times within this window (condition (JTJ) ensures that no input impulse is lost). BN 
receives impulses both from the Poisson stream and from the line. But no more than 
one impulse from the line may have time to reach BN's input during time interval less 
than A. Therefore, the other (2n + 1) impulses must be received from the Poisson 
stream. On the other hand, if 2(n + 1) input impulses are received from the Poisson 
stream during time interval ]0; A[, the inequality ( fl8l) holds for sure, no matter is the 
impulse from the feedback line involved, or not. Therefore, pa > p{2n + 2, A) > 0, 
where p(i, A) gives the probability to obtain i impulses from the Poisson stream during 
time interval A [16J: p(i, A) = e~ AA (AA)7^!. 

Having in mind f fl8|) . let us split the integration domain for sq in (jSJ) in the following 

way: 

/ ds = ds + > / ds + / ds , 

and introduce the notations: 

E5= % a a n+1 



h= / ds dsi... / ds n+ iP(t ,s ) |"^[ P(t k , s k \ t k -i, s k -i), 

z = 0,l,2,. ..,n, (19) 
/ ds / dsi . . . / ds n+ xP(t ,s ) Y[ p (tk,s k \t k -i,s k -i), (20) 



3=0 



fe=l 
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where we assume, that YljLji = ® f° r h > h- 

Consider the fixed sequence of ISIs, (t , . . . , t n ), which belongs to D x . Domain of s 
values covered by Ij, i = 0, . . . , n, corresponds to the scenario, when impulse, which was 
in the feedback line at the beginning of interval to (with time to live Sq), will reach BN 
during interval £j. In this process, after each firing, which starts ISI t k , k < i, the time 
to live of the impulse in the feedback line is decreased exactly by tk-i- This means, that 
variables of integration {s , . . . , s n+ i}, above, are not actually independent, but must 
satisfy the following relations: 

k-i 

Sk = s - 1 ^2t j , k = l,...,i, (21) 

j=0 

which are ensured by 5-function in the bottom line of ([IT]) . Next to time to live must 
be equal A: 

s i+ i = A, (22) 

and this is ensured by ^-function in the top line of (fT7l) . The next to Sj + i times to live 
again are decreased by corresponding ISI with each triggering. Due to (1T8]) . this brings 
about another set of relations: 

k-l 

s k = A-J2 t h k = i + 2,...,n + l, (23) 
j=i+i 

which are again ensured by 5-function in the bottom line of (IT7I) . Relations Q21j) . (T22l) 
and ( 123|) together with limits of integration over sq in ( fl9|) ensure that at D\ the following 
inequalities hold: 

s k >t k , fc = 0, 

Si < U, (24) 
s k >t k , k = i + l,...,n. 

Inequalities fT24"j) allow one to decide correctly which part of rhs of ffTTj) should replace 
each transition probability P(t k ,s k \ t k -i,s k -\) in ( fl9|) . and perform all but one 
integration. This gives: 



h= / ds / dsi ■ . . . ■ / ds„+iF(t | s )f(s ) F(t fc | s k )5(s k - s + ^ t 

J J J i — i ■ — n 



k-l 

'J J 



x^i-if. 



k=l 3=0 
n+1 Jfe-1 

x I s i+1 ) 8{s i+1 - A) JJ | s k )5(s k - A + ^ tj) 

k=i+2 j=i+l 
n n—l 

= F(t n+1 | A - J2 tiWn I A - *i) ' • • • • F (^+ 2 I A - U+i)F(t l+1 | A) 

.7=i+l i=i+l 
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E y* J i-l i-2 

x / F{U | s - ^2 tj)F{U-x \ so-^tj) ■ F(t x | s - t )F(t | s )/(s )ds , 

^ i . i=o 3=0 

2^=0 1 j 

i = 0,1,2,..., n. (25) 

The last expression might be obtained as well by means of consecutive substitution 
of either top, or bottom line of f|T7|) into f|T9|) . without previously discovering f l2Tj) - 

(El. 

Finally, integral I n+ i corresponds to the case, when at the beginning of interval 
t n+ i, the line still keeps the same impulse as at the beginning of t Q . Therefore, I n+ \ 
comprises the rest of scenarios contributing to the value of P(t n+ i, . . . , t ) in (Q. Here, 
the bottom line of (jTTJ) ensures that values of variables of integration {so, . . . , s n+ i}, 
which contribute to the I n +i, should satisfy the following relations: 

fc-i 

Sk = s -^tj, k = l,...,n + l, (26) 

3=0 

which taken at the domain Di, defined in ( ITBl) . results in inequalities 

s k >t k , k = 0,...,n. (27) 
Equations (ITT)) and fT27|) allow one to perform integration in (120]) and to obtain: 

n+l 



4+1 = / ds / dsi . . . / ds„ + iF(t | s )f(s ) TT F(t k \ s k )5{s k - s + V) t 
^E^o^ -/o fc=1 i=0 ' 

/n n— 1 

^(Wi I so - i s ° - Yl ^ ■ ■ ■ F ( fi i s ° - to ) 

v t 3=0 3=0 

x F(t I s )/(so)ds . (28) 

Taking into account (125)) and f )28|) . one obtains the following expression for joint 
probability density P(t n+ i, . . . ,t ): 

n+l 

P(tn+1, ■ ■ ■ , to) = h 
i=0 

n n+l fc— 1 

=^F(t m |A) n F(t fc |A- 

8=0 fc=j+2 j=i+l 

J"" ' ^0 | S ) /(S ) I]X^ I S ° - J2^ dS0 

^)=o *i fc=l j=0 

/A n+l ft— 1 n 

F(* | so) /(so) H I s ° - X^) ds °> $^^<A,(29) 

^j=o*i fc=l j=0 i=0 

where we assume, that J^L^ = an d HiL^ = 1 f° r Ji > h- 
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Using fl7|), for conditional probability density P(t n+ \ \t n , . . . ,t ) one derives: 

n+l fe-1 

X 



P(t n+l | t n , . . . ,to) = 1 • ($_ F(* <+1 | A) J] F(t k \A-Y, l i 

r\i n ,...,i ) v. =Q k=i+2 j= . +1 

/v*._ t- * i 

^}=o'j fc=l j=0 

/A n+l fc— 1 

F(t | s ) /(s ) J] F(t k \s -J2 h) d 

^™=0*i fc=l i=0 



„v^ 1 / « fc— 1 

X 

Sj=o *J k=l j=0 

„/\ n+l fc— 1 

' ^ l.s,, 

j=0 

n 

J2U<A, (30) 

i=0 

where expression for P(t n , . . . ,t ) can be obtained from (1291) with {n — 1) substituted 
instead of n. 

5-4- Singular part of P(t n+1 , . . . , to) 

In order to obtain the singular part of expression, defined in (129]) . let us first derive 
singular parts for all Jj, i — 0, . . . , n and I n+1 separately. In order to keep the expressions 
shorter, we represent Ii as follows 

Ii(to, ■ ■ ■ , t n +\) = Xi(to, . . . ,ti) ■ Yi(ti + i, . . . , t n +i), i = 0, 1, . . . , n, (31) 

where 



i-l i-2 



F(ti\s -J2tj) F (ti-i\ s a ~ ■ -F(ti\sa - t )F(t \s )f(s )ds , (32) 

— n A — n 



y>i-l 4 . 



j=0 j=0 



n-1 



F< = F(t n+1 \A- J2 hWitn I A - 53 tj) . . . F{t i+2 | A - t i+l )F{t l+1 | A). (33) 

j=i+l j=i+l 

It is clear, that at D\, Xi is the part of the probability density for (i + 1) successive 
ISI durations, which corresponds to the case when the impulse, which was in the line 
at the beginning of the first ISI, reaches the neuron's input within the last one. And 
the Yi gives the probability density for (n + l —i) successive ISI durations provided the 
impulse enters the line just at the beginning of the first one of these ISIs. 

At the domain considered, namely, for Yl7=o^ < ^ the expressions for F(t n \ 
A — YujZi+i tj)> • • • > F(t i+2 | A — t i+ i) and F(t i+ i | A) have no singularities, see ( flOl) . 
Therefore 

n n— 1 

F sing = F smg (Wi | A _ ^ tj )F(tn\±- £ tj) . . . F(t i+2 \A - t i+1 )F(t i+l \A). (34) 

j=i+i j=i+i 

_____ siii£C 
At the same time, intergation limits in (1321) ensure that X i 6 = 0. Indeed, each 

integral Xi (and, originally, !»), i = 0,1,..., n, covers the half-open interval so G 



X 
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] X/j-=o*J' X/j=o^']' ^ ne on ^ singularity of integrand in (I3"2"j) at this domain is 
5(X)}=q*j — s o) provided by | s — 5^i=o*i)> see P^ - an d it disappears after 

intergation. Therefore 

n n—1 

J sing = FS ing (Wi|A _ £ t])F{tn \ A _ *.)...F(t i+2 |A-t m )F(* i+1 |A) 

j=i+l j=i+l 
i- i— 1 i— 2 

f ~" ' Ffe I so -^t^Ffc-i I so - • I s ° - t )F(t I s )/(so)ds , 

^}=o t i j=0 j=0 

i = 0,l,...,n. (35) 
In the same way, for the singular part of integral I n+ i one obtains 

n n—1 

J^ g = a • F sin g(t n+1 | A - ^)F(* n | A - ^) . . . Ffa | A - t )F(t Q | A), (36) 

j=0 j=0 

where a is the (^-function's mass in f(s), see ( TTTT) . 

Taking into account (ITDj). ( 135]) and f l3"6"j) . for the singular part of the probability 
density P(t n+ i, . . . , t ) one obtains 

n+l 

P sin g(t n+1 ,t n ,...,t ) = ^/r g 

i=0 

n n+l 



A i ■ <*( X *J - A ) + • + • • • + t - A), 



i=0 j=i+l 

n 

J><A, (37) 

8=0 

where A4 and A n+1 denote regular factors, defined by the following expressions: 

n-l 

Ai = Xt n+1 e- At " +1 ■ F(t n I A - ^t 3 )... F(t i+2 | A - t l+1 )F(t i+1 | A) 

3=i+l 

/y>*_ t- i— 1 i— 2 

3 I s o -^y^i-i I s - • I s o ~ I s o)/(s )ds , 

^j=o *j j=0 j=0 

i = 0,l,...,n, (38) 

n-l 

A n+1 = a • Xt n+1 e~ xt ^ • F(t n I A - X t,) . . . F(h | A - * )F(* I A). (39) 

3=0 

Obviously, each factor Ai, i = 0, ...,n, gives the probability to obtain (n + 1 — i) 
successive output ISIs of overall duration exactly equal A. And A n+ i gives the 
probability to obtain (n + 2) successive output ISIs of overall duration exactly equal A. 

The presence of 5-functions in joint probability density P(t n+ i, . . . , i ) can be 
additionally explained as follows. If at the beginning of (i + l)-th ISI, the impulse enters 
the line, then output interval t n+ i will start with that same impulse in the feedback line 
with time to live equal s n+ \ = A — 52™_ i+1 t/. To trigger BN after time exactly equal 
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s n+ i after that, it is nesessary to obtain one impulse from the Poisson stream during 
time interval s n+ i. This event has non-zero probability, therefore we have the non-zero 
probability of an output ISI exactly equal to s n+ ±: t n+ i = A — Y^j=i+i tj- This gives the 
corresponding (^-functions in ISI probability density. The term with S(t n+ i + . . .+t — A) 
corresponds to the case, when the impulse enters the line at the beginning of to- 

From (J7|) and f )37|) one can easily derive the following expression for the conditional 
probability density: 

P-S( f „ +1 K,...,to)= ■ ¥K ^— ) £> .,(£«,- a) 



j=i+l 



A n 
+ p( . " +1 . s -S(t n+1 + ... + t -A), $><A, (40) 
f[t n ,...,t Q ) 

where A4 and A n+ i are defined in f l38jl and fl39l) . It should be outlined, that joint 
probability density P(t n , . . . , to) has no singularities at the domain t n < A — Y^h=o 
see (137|) with (n — 1) substituted instead of n. 

As one can see, function P(t n+1 \ t n , . . . ,t ) contains singularty at t n+ \ = A — t n — 
t n -\ — ... — t - The dependence of the singular part of function P(t n +i \ t n , . . . , t ) on 
to cannot be compensated by any regular summands, therefore, the whole conditional 
probability density P(t n +i | t n ,...,t ) depends on to- It means, that the condition 
(|2J) does not hold for any n for the output stream of BN with delayed feedback. The 
Theorem 1 is proven. 



6. Particular cases 



In previous section, we have prooven the impossibility to represent the stream of output 
ISI durations for BN with delayed feedback as a Markov chain of any finite order. In 
particular, output ISI stream is neither a sequence of independent random variables, 
and therefore is non-renewal, nor it is the first-order Markovian process. 

In the course of proving Theorem 1 (see Sections H] and [SJ, we have obtained the 
expression for P(t n+ \ \ t n ,...,t ) at the domain Y^i=o^i < A in general case of an 
arbitrary n, see (1301) . 

In this section, we consider the two particular cases of P(t n +i | t n ,...,t ) when 
n = and n — 1, namely, the single-moment conditional probability density P(ti | to) 
and the double-moment conditional probability density P(t 2 | ti,to) and obtain the 
expressions for P(ti | t ) and P(t 2 | ti,t ) for domain ( fl8i) . as well as for all other 
possible domains, which were omitted in general consideration. 



6.1. Conditional probability density P(t x | t ) 

In order to derive the exact expression for conditional probability density P(ti | t ) for 
neighbouring ISI durations, we take Steps 1-3, outlined in Section HI for n = 0. In the 
case of P(ti I to), there are only three domains, on which the expressions should be 
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obtained separately, namely cases to < A, t > A and t = A. Performing intergation 
in OH]), one obtains the following expressions for P(ti,t ) at these domains: 

P(t u to) = F{tt | A)P(to), t > A, (41) 

= F{h | A) [° F(t | so)/(s )d So 

+ / F{U | s - to)F{to | so)/(s )ds , £ < A. (42) 



./to 

Then, by definition of conditional probability densities, one obtains: 

P{h | t ) = F{t x | A), to > A, (43) 

= — { F (ti I A) J ° F(t | s )f(so)d So 

+ / F{tt | s - t )F(t | So)f(s )dso) , t < A. (44) 

Jt ' 

It should be outlined, that the output ISI probability density P(t Q ) has no singularities 
at the domain t < A. Indeed, due to f lT0|) -f lT2|) . the only 5-function contained in -P(to) 
is placed at t = A, see Figure [3] (a). 

In vicinity of the point t = A, the single-moment conditional probability density 
can be derived as 

/ dt P(ti,t ) 

Pit, | t = A) = [im^ , (45) 

/ dtoP(to) 

A-e 

which just gives 5-functions' masses both in numerator and denominator, and delivers 
Pit, | t ) = F(ti | A), t = A. (46) 

Expressions (|43|) . ( jH| and ( 146|) can be understood as follows. Since t > A, one 
can be sure that the line has time to get free from impulse during to, therefore at the 
moment of next firing (at the beginning of t\) the impulse enters the line and has time 
to live equal A. In the case of to < A, see (jSJ), two possibilities arise. The first term 
corresponds to the scenario, when the feedback line discharges conveyed impulse within 
time interval to, and the second one represents the case when at the beginning of t\ the 
line still keeps the same impulse as at the beginning of to- 

It can be shown, that the following normalization conditions take place: 

oo oo 

/ dtiP(tx | t ) = 1, and J dt P(t!,to) = P(ti). 



The singular part of P(ti | t ) can be easily extracted: 

P sing {t 1 \t ) = e- XA XA-5{t l -A), t >A, (47) 

Ati e~ Atl ' rt0 



- ( j ° F(t | s )f(so)dso ■ 5(h - A) + 
+ a F(t | A) • 5(t + ti - A)), t < A. (4? 
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Figure 4. Conditional probability density P(t\ \ to) for r = 10 ms, A = 8 ms, A 
= 150 s _1 , iVo = 2, to=6 ms (a) and to= 11 ms (b), found numerically by means of 
Monte-Carlo method (the number of firings accounted N = 30 000). 



Obviously, expression (HHl) could be obtained directly from (|38l) -(H0l) by substituting 

71= 0. 

As it can be seen from (H71) and ( 1481) . the number of (^-functions in P{t\ \ to) and 
their positions depend on t , therefore the conditional probability density P{ti \ t Q ) 
cannot be reduced to output ISI probability density -P(ti). Therefore, the neihgbouring 
output ISIs of BN with delayed feedback are correlated, as expected. 

Examples of P(ti \ t Q ), found for two domains numerically, by means of Monte-Carlo 
method (see Section [7] for details), are placed at Figure HI 



6.2. Conditional probability density P(t 2 \ ti,t ) 

In order to derive the exact expression for conditional probability density P{t 2 \ £i,to) 
for the succecive ISI durations, we take Steps 1-3, outlined in Section HJ for n = 1. 
In the case of P(t 2 ,ti,to), there are six domains, on which the expressions should be 
obtained separately, namely, the domain 

D x = {tx,t | h + t < A}, 

which was already utilized in Section [51 and five remaining: 



D 2 


= {h 


to I 


to 


> A 


and 


ti 


> A}, 


D 3 


= {ti 


to | 


to 


< A 


and 


ti 


> A}, 


D 4 


= {h 


to I 


to 


> A 


and 


h 


< A}, 


D 5 


= {ti 


to 1 


to 


< A 


and 


A 


-t < 


d 


= {h 


to 1 


to 


+ h = 


A}. 







In the case, when the exact equality to + ^i = A holds, namely, if (ti, t ) G d, the product 
P(t 2 | ti,to)dt 2 gives the probability to obtain an output ISI of duration within interval 
[t 2 ;t 2 + dt 2 [, provided the overall duration of two previous ISIs accurately equals A. 
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Expressions for P(t2 \ ti,t ) can be found exactly on each domain: 

P(t 2 \t 1 ,t ) = F(t 2 \A), (Mi)eAj, (49) 

= F(t 2 \A), (t ,t x )eD 3 , (50) 

= F(t 2 \A), (t ,h)ed, (51) 

= F{t 2 \A-t 1 ), (to,ti)eD A , (52) 

= p ^ to) (P(*2 | A - tjFfa | A) J ° F{t | so)/ (so) ds 

+ F(t 2 \A) [ F(t 1 \s -t )F(t \s )f(s )ds ), (t ,h)eD 5 , (53) 

./to 7 

= P // I \ f^(*2 I A - tjFih | A) / ° F(t | s )/(s ) ds 

-r(.Cl ) t 0j V JO 

+ F(t 2 | A) f F(tx | s - t )F{t | so)/(s )ds 

J to 

+ / F{t 2 \s -t -t 1 )F{t 1 \s -t )F{t \s )f{s )ds ), 

{t 0i t 1 )eD 1 . (54) 

where P{h, t ) = F(t x \ A) f*° F(t \ s )f{s )ds + f*F{t 1 \ s -t )F(t | s )f(s )ds , 
according to (jUj). 

The probability density P(ti,t ) contains <5-function at the domain d, see (|4~BI) . In 
(|5Tj) . the two-time conditional probability density was derived as 

A-to+e 

j dfiP(t 2 ,*i,to) 
-Pfe I ti,t ) = 1 ™ A ~ f °~ £ . , (*o,*i)ed, 

e— >0 ii— co+e 

/ dtiP(ti,t ) 

A-to-e 

compare with (1451) . 

It is worth to notice, that P(ti,t ) is regular function on both Di and D$, see (1531 
and fl54l) . Indeed, from ({IT]) and (j4"81) one can see, that P(ti, to) ma Y include singularities 
only at the points ti = A and t\ — A — to- None of these points belongs to Di, or D 5 . 

It can be shown, that the following normalization conditions take place: 

oo oo 

/ dt 2 P(t 2 I t 1; t ) = 1, and / dt P(t 2 ,t u t ) = P(h,tx). 



The singular part of the conditional probability density P(t 2 | ti, to) can be derived 
as follows: 

P sing (t 2 | h, t ) = e~ A ' 2 At 2 • 5{t 2 - A), (t , t x ) ED 2 UD 3 Ud, (55) 

= e" A ' 2 At 2 • + t 2 - A), (t , t x ) e D A . (56) 

= WTTT • ( F ^ I A ) / F ^ I s ^f^ ds ° • 6 ^ + *2 - A) 

-r^lj^oj v Jo 

+ / F(ti | s - t )P(t | so)/(s )ds • 5(t 2 - A) 
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Figure 5. Conditional probability density P(ta | ii,io) f° r t = 10 ms, A = 8 ms, A 
= 150 s _1 , Nq — 2, ti=13 ms, trj=13 ms (a) and ti = 6 ms, to = 13 ms (b), found 
numerically by means of Monte-Carlo method (N — 30 000). 



to+tl 



(t ,*i) e As, 



(57) 



^(ti s *o) 

+ | A) I ° F(t | so)/(s )ds • 5(h + 1 2 - A) 
./o 

+ a ■ F(t! | A - t )F(t | A) • S(t + t 1 +t 2 - A)) , 

(<o,*i)e£>i. (58) 

Obviously, expression (1581) could be obtained directly from (|38l) -( l40l) by substituting 
n — 1. 

As one can see, the singular part of P(t 2 | ti,t ) depends on to, therefore P(t 2 | ti, to) 
cannot be reduced to P(t 2 \ ti), which means that the output stream is not first-order 
Markovian. 

Examples of P(t 2 \ ti,to), found numerically for different domains, are placed at 
Figures [5] and [61 



7. Numerical simulation 



In order to check the correctness of obtained analytical expressions, and also to 
investigate wheather the output ISIs stream is non-Markovian for BN with higher 
thresholds as well as for N = 2, numerical simulations were performed. A C++ 
program, containing class, which models the operation manner of BN with delayed 
feedback, was developed. Object of this class receives the sequence of pseudorandom 
numbers with Poisson probability density to its input. The required sequences were 
generated by means of utilities from the GNU Scientific Librar}|§| with the Mersenne 
Twister generator as source of pseudorandom numbers. 

§ http : / /www . gnu. org / software /gsl/ 
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Figure 6. Conditional probability density P{t% \ ti,to) for r = 10 ms, A = 8 ms, 
A = 150 s _1 , Nq = 2, ti—Q ms, to— 3 ms (a) and t\ = 6 ms, to = 1 ms (b), found 
numerically by means of Monte-Carlo method (N — 30 000). 



Program contains function, the time engine, which brings system to the moment 
just before the next input signal, bypassing moments, when neither external Poisson 
impulse, nor impulse from the feedback line comes. So, only the essential discret events 
are accounted. It allows one to make exact calculations faster as compared to the 
algorithm where time advances gradually by adding small timesteps. 

The conditional probability densities, P(t 1 \ t ) and P(t 2 \ t 1 ,t ), are found 
by counting the number of output ISI of different durations and normalization (see 
Figures H] - Ej). Obviously, for calculation of conditional distiributions only those 
ISIs are selected, which follow one or two ISIs of fixed duration, i for P{t\ \ t ) 
and {t\,to} for Pfa | £i,to)- The quantity, the position and the mass of delta- 
functions, obtained in numerical experiments for BN with threshold 2, coincide with 
those predicted analitycally in fl47|) . f|48j) and ([55]) - ( |58l) . 

For iVo > 2, conditional probability densities P{t\ \ to) and P{t 2 \ ti,to) are similar 
to those, found for N =2. In particular, both the quantity and position of delta-functions 
coincide with those obtained for BN with threshold 2, as expected, compare Figures [7] 
and |6j 

8. Conclusions and discussion 

Our results reveal the influence of the delayed feedback presence on the neuronal 
firing statistics. In contrast to the cases of BN without feedback [17] and BN with 
instantaneous feedback [18], the nighbouring output ISIs of BN with delayed feedback 
are mutually correlated. It means that even in the simplest possible recurrent network 
the ISI stream cannot be treated as the renewal one. The presence of nearest ISIs 
correlation was reported for spike trains of a neurons in different CNS and peripheral 
NS structures [T9l [20]. 

Moreover, we prove, that the output ISI stream of BN with delayed feedback cannot 
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Figure 7. Conditional probability density Pfa | ii,io) for r = 10 ms, A = 8 ms, 
A = 800 s -1 , A^o = 4, ti—6 ms, io=3 ms (a) and t\ = 6 ms, to = 1 ms (b), found 
numerically by means of Monte-Carlo method (N — 30 000). 

be represented as the Markov chain of any finite order. This is in accordance with rare 
attempts of experimental estimation of the Markov order of neuronal spike trains (see, 
e.g. [3], where it is established that the order, if any, must be greater than 3). 

We expect the same non-markovian property for firing statistics of any single neuron 
with delayed feedback, whatever neuronal model is used, and conclude that it is namely 
the delayed feedback presence results in non-markovian statistics found. One should 
take this fact into account during analysis of neuronal spike trains obtained from any 
recurent network. 
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